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ABSTRACT 

We have carried out a numerical study of the effect of large scale magnetic 
fields on the rate of accretion from a uniform, isothermal gas onto a resistive, 
stationary point mass. Only mass, not magnetic flux, accretes onto the point 
mass. The simulations for this study avoid complications arising from bound- 
ary conditions by keeping the boundaries far from the accreting object. Our 
simulations leverage adaptive refinement methodology to attain high spatial fi- 
delity close to the accreting object. Our results are particularly relevant to the 
problem of star formation from a magnetized molecular cloud in which thermal 
energy is radiated away on time scales much shorter than the dynamical time 
scale. Contrary to the adiabatic case, our simulations show convergence toward 
a finite accretion rate in the limit in which the radius of the accreting object 
vanishes, regardless of magnetic field strength. For very weak magnetic fields, 
the accretion rate first approaches the Bondi value and then drops by a factor 
~ 2 as magnetic flux builds up near the point mass. For strong magnetic fields, 
the steady-state accretion rate is reduced by a factor ~ 0.2/3 1 / 2 compared to the 
Bondi value, where (3 is the ratio of the gas pressure to the magnetic pressure. 
We give a simple expression for the accretion rate as a function of the magnetic 
field strength. Approximate analytic results are given in the Appendixes for both 
time-dependent accretion in the limit of weak magnetic fields and steady-state 
accretion for the case of strong magnetic fields. 

Subject headings: ISM: magnetic fields — magnetohydrodynamics (MHD) - 
stars: formation 
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Introduction 



Accretion of a background gas onto a central gravitating body is of central importance in 
astrophysics. Examples range from protostellar accretion from molecular cores to accretion 
of interstellar gas in galactic nuclei. The classical late-time solution for the case of a central 
point of mass M* immersed in a uniform, initially static, unmagnetized gas was given by 



Bondi (1952) as 



4n\r B 2 p c 
GM* 



(1) 
(2) 



where and are the sound speed and density of the background gas, Mb is the steady- 
state rate of accretion onto the central particle, tb is the Bondi length which characterizes 
the dynamical length of the inflow and A is a dimensionless parameter that depends on the 
equation of state of the background gas. For the isothermal case, A = exp(1.5)/4. The 
Bondi time t B = tb/cqo defines the dynamical time for this accretion process. This basic 
model has been extended to more general cases by numerous authors. These generalizations 



include non-stationary central particles (Bondi & Hoyle 1944 Shima et al. 1985 Ruffert 



1994 Ruffert & Arnett 1994), the cases of ambient gas with net vorticity (Krumholz et al. 



2005), turbulent ambient gas (Krumholz et al. 2006), magnetized accretion from ambient 



gas threaded by both large (Igumenshchev & Narayan 2002 Pang et al. 2011) and small 



(Shapiro 1973 Igumenshchev 2006) scale magnetic field topologies, the case of a turbulent, 



magnetized ambient gas ( Shcherbakov||2008 ) , and the case of accretion onto magnetized stars 



(Toropin et al. 1999 Ustyugova et al. 2006 Kulkarni & Romanova 2008 Romanova et al. 



2008 Long et al. 2011; Romanova et al. 2011), to name a few. 



Stars form via gravitational collapse, at least initially ( McKee &: Ostrik"er||2007 ). There- 
after, gas may accrete onto the star from the ambient medium. If the star has a supersonic 



motion relative to the ambient medium, this subsequent accretion is negligible (Krumholz et 



al. 2005), but if the star is moving slowly, the accretion can be significant, which forms the 



basis for the competitive accretion model for star formation (e.g., Bonnell et al.|199~ ). There 
exists ample evidence that the gas in molecular clouds and cores is threaded by strong mag- 



netic fields (Crutcher 1999 McKee & Ostriker 2007). Furthermore, star forming molecular 



clouds are well characterized as "radiatively efficient" in that gas heating due to compres- 
sional motion is rapidly radiated by thermally excited dust and molecules. These consider- 
ations thus motivate the study of Bondi-type accretion of an isothermal gas threaded by an 
initially uniform magnetic field onto a point mass. We address this problem with the RAMSES 



magnetohydrodynamic (MHD) code (Teyssier 2002) and conduct a parameter study over a 



range of magnetic field strengths thought to be relevant to star formation. Our simulations 
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leverage the adaptive mesh refinement (AMR) capability of the code to retain high spatial 
resolution close to the accreting object while keeping the boundaries of the computational 
domain far from the accreting object. We discuss the results of mesh convergence studies 
and compare our numerical results against analytic calculations in the limiting case of a 
dynamically weak magnetic field to verify our calculations. We also compare our numerical 
results against simple analytic approximations in the case of a strong magnetic field. 



2. Numerical Setup 

Our numerical models consist of a Cartesian computational domain that extends from 
— 25re to 25re in each direction. The domain is initialized with an isothermal, perfectly 
conducting, uniform collisional gas with initial magnetic field in the z direction. We consider 
the cases with an initial thermal to magnetic pressure ratio, = 8vrP /B 2 , of 1000, 100, 10, 
1, 0.1 and 0.01. The RAMSES code has been used to evolve this state forward according to 
the equations of ideal, isothermal MHD, 

^ + V-pv = -S M (3) 
9 " V ■V.(pvv) + vfp + g!V (B - V)B = -g*g* (4) 



dt vr ' \ 8n J 4vr x 2 

^--Vx(vxB) = (5) 

P = pc\ (6) 

where p is the gas density, v is the velocity, B is the magnetic field, P is the thermal pressure 
and c is the isothermal sound speed. These equations include the gravitational force due to 
a point particle of mass M* of F g = —GM*p x/x 2 . 

The key assumption we make in our treatment is that the point mass accretes mass, 
but not flux. Observations show that the magnetic flux in young stellar objects is orders of 
magnitude less than that in the gas that formed these objects, implying that flux accretion 



is very inefficient, presumably due to non-ideal MHD effects, including reconnection (McKee 
fc Ostriker|[2007| . We model mass accretion onto the central point mass by including a mass 
sink term but no flux sink term inside a radius, r acc = 4Aa;, equal to four grid zones on the 
finest AMR level. The effect of the accreting particle is coupled to the dynamical equations 
through the source term, 

i -J-max(p-— ? ,0] if |x| < r acc 

3 otherwise, 
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where At is the time step on the finest AMR grid level and v A,max is the maximum Alfven 
speed, I?/(47rp) 1//2 , within a radius of 6Ax around the accreting particle. Under this con- 
struction, the accreting particle absorbs all but enough of the mass entering the accreting 
particle radius so that the local Alfven speed never exceeds tA.max- Thus, the accreting 
particle always absorbs the largest quantity of mass in the local region possible without 
introducing new local extrema in the Alfven speed. This prevents the accretion source from 
imposing a stringent (or vanishingly small) constraint on the maximum numerically stable 
time-step at the expense of some artificial clipping of the Alfven speed in the inner few zones 
around the accreting particle. We note that in all of the models considered in this paper, the 
initial gas density is sufficiently low that the total mass accreted onto the central particle is 
negligible compared to M* and that self-gravity in the ambient medium may be neglected. 

We discretize the numerical domain onto a base level grid of 64 3 . For the purposes 
of describing the initial mesh we will denote this level as I = 0. We note, however, that 
the RAMSES AMR implementation uses an oct-tree data structure for level traversals that 
always denotes level indices by the base 2 logarithm of their resolution. In our models, 
/ramses — l°g264 + 1 = 6 + 1. Successive levels are chosen for refinement by an increment of 
2 3 in grid zone density in a geometrically nested fashion according to the criterion 

25r B 

n < -2T, ( g ) 

where r\ indicates the radius of the spherical refined region on the level I. We further impose 
the additional criterion that any zones containing steep density gradients Vp • Ax/p > 1/2 
are also refined, independent of location. This second refinement criterion is met only at 
late times after non-axisymmetric flow patterns have set in, and it triggers only on transient 
flow features. Most of the models were refined to a maximum level I = 8 for an effective 
resolution of 64 x 2 8 /50 ~ 328 zones per thermal Bondi radius on the finest level. In the 
cases with a strong initial magnetic field, it is also useful to consider the numerical resolution 
on the scale of the 'Alfven-Bondi'' radius, 

GM * 1 R fQ\ 

v A z 

The finest mesh resolution per Bondi radius, mesh resolution per Alfven-Bondi radius and 
total simulated time for each model is tabulated in Table [Tj We note that the magnetic length 
scales are well resolved for all but the case of /3 = 0.01. We therefore will consider only the 
models with /3 > 0.1 for the majority of the analysis presented in this paper. The (3 = 0.01 
model is used only to extract an estimate of the steady accretion rate over a wider range 
of magnetic field strengths. We do note, however, that numerical mesh convergence studies 
have shown our models to be within the range of asymptotic convergence with a Richardson- 
extrapolation error estimate on the average accretion rate of 14% or less at late times. A 
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detailed discussion of the numerical convergence properties of our models is presented in 
Appendix [Aj Each of the models were run to a final time t e nd sufficiently long to attain a 
statistically steady accretion rate onto the central particle. 

3. Results 
3.1. Morphology 

We begin by discussing the gross morphological flow features and their development for 
each of the numerical models. These flows are well illustrated by slices in the y-z plane of 
density, the direction of magnetic flux and velocity as shown at several times for each model 
in Figure [TJ Initially parallel magnetic fields are amplified as they are dragged inward by 
the global accretion flow, eventually suppressing accretion in the equatorial plane. Inflow 
along magnetic field lines, on the other hand, is uninhibited by magnetic pressure. This flow 
configuration leads to the evacuation of gas from the poleward directions into a thin, dense, 
irrotational disk in the midplane. 

Accumulation of mass in the midplane is accompanied by a corresponding increase in the 
inward gravitational attraction. The magnetic flux tubes that thread the disk are gradually 
pulled further toward the accreting particle as the accumulation of mass in the midplane 
continues. We support this picture more quantitatively in Figure [2] We use w to denote the 
cylindrical radius and plot the ratio of the mass influx in the equatorial direction 

$m,vj = / pv • w sin 8d6d(j) (10) 
Js 



Table 1. Simulation Parameters. 
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Fig. 1. — Slices in the y-z plane showing the inner (2r^) 2 of the numerical models with initial 
magnetic field strengths of /3 = 100, — 10, /3 — 1, and /3 = 0.1 from top to bottom. The 
state of the numerical models are shown at the times t = 0.5, 1.5, and 3.0 £b from left to 
right. The colormap indicates log 10 (p/p o ), green lines represent magnetic flux tubes drawn 
from equidistant foot-points in the midplane and the white arrows indicate the flow pattern 
in the plane of the slice. 
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Fig. 2. — The ratio of the mass influx in the equatorial direction to the mass influx in the 
polar direction for several magnetized models at t — tend; scaled so that uniform spherical 
inflow takes a value of unity. 

to the mass influx in the polar direction 

®Mz = I P v • z sin OdOdcf) (11) 



Js 

along a spherical control surface S of radius r for each of the magnetized models at t — t e nd- 
The curves in Figure [2] have been scaled by a constant 2/n so that uniform spherical inflow 
takes a value of unity At large distances (|x| > re), the flows become increasingly dominated 
by polar inflow with increasing initial magnetic field strength. However, at smaller distances 
(|x| < re), the cylindrical to polar influx asymptotes toward ~ 2 with increasing magnetic 
field strength. On smaller scales where magnetic forces break spherical symmetry, the mass 
influx is predominantly along the equator. 

As infall in the midplane proceeds, flux tubes that reach the accreting particle are instan- 
taneously liberated from the accreted mass and accompanying gravitational force anchoring 
them. This causes episodic releases of strong, outward propagating flow. This configuration 
of outflow driven by magnetic buoyancy is known as the magnetic interchange instability 



(Bernstein et al. 1958; Furth et al. 1963). In the models with moderate or strong initial 
magnetic fields strengths, corresponding to (3 = 10, (3 = 1 and (3 = 0.1, interchange unstable 
flows originating at r acc <C t-q lead to episodes of net outflow out to radii comparable t-q 
in the equatorial plane. Flux tubes that are outwardly released by resistive accretion are 
prevented from escaping completely by the continued accretion pressure of the surrounding 
gas. The net mass inflow in these models is therefore mediated by the rate at which inflow- 
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ing material percolates through this non-axisymmetric network of magnetically buoyant flow 
close to the accreting particle. 

The models attain magnetic forces that balance F g at r ~ re/2 in the midplane by the 
time steady accretion sets in, independent of the initial (3. The weak magnetic field lines 
in the = 100 case become highly stretched before they are strong enough to provide any 
resistance to being swept further inward as shown in the top row of Figure [TJ This flow 
leads to the development of strong, thin current sheets and oppositely directed magnetic 
field lines that closely approach each other in the midplane. This configuration is unstable 



to reconnection in magnetic resistive tearing modes (Furth et al. 1963 Rutherford 1973). 
In the case of our numerical code (and all ideal MHD codes), resistive reconnection occurs 
when oppositely directed magnetic flux tubes become separated by < Ax and unresolved. 
While the size scale of the "magnetic islands" generated through this process is determined 
by the numerical zone size, our numerical resolution is adequate to be sure that this size 
scale is small compared to dynamical scale associated with thermal (Ax <C tb) and magnetic 
(Ax <C tab) force gradients. Furthermore, we have carried out resolution studies to ensure 
that the resoultion used in our models is sufficient to yield a converged late-time accretion 
rate. Ultimately, mass inflow is limited by the rate of production of magnetically isolated 
islands by tearing mode reconnection in regions characterized by thin, strong current sheets. 
These islands continue toward the accreting particle, unconnected to the global magnetic 
field structure. As a means to visualize flows that are most susceptible to reconnection by 
numerical resistivity, we define the magnetic shear parameter 

= Ax-(VxB) 

Amag i-gi • l xz V 

Regions near or exceeding a magnetic shear parameter of ~ 1 are highly susceptible to 
reconnection via magnetic tearing modes. Figure [3] gives a three dimensional sense of the 
geometry and scale of the flows subject to numerical reconnection by plotting isosurfaces 
of the magnetic shear parameter at t — t cn d, indicating efficient numerical reconnection on 
scales of r < re/2. Reconnection events release magnetic tension that leads to magnetically 
tangled, non-axisymmetric flow in this region. 



3.2. Comparison to Analytic Predictions for High f3 Flow 

Analytic predictions of the behavior of the accretion flows for the limiting case of dy- 
namically weak magnetic field are derived in appendix [B} The focus of this section is to 
compare the results of the /3 = 100 numerical model with these analytic predictions. Equa- 



tions (B16), (B23) and (B25) give predictions of the steady state gas density, radial magnetic 



Fig. 3. — Isosurfaces showing the innmermost (1.5rB) 3 of the magnetic shear parameter Xmag 
at t — t eri( i for the f3 = 100 model indicating regions of magnetic reconnection due to tearing 
mode instability. Blue curves represent magnetic field lines with footpoints evenly spaced 
along the y coordinate axis. 



field and non-radial magnetic field respectively. (Results for the accretion rate will be dis- 



cussed in £3.3) It should be emphasized that r in these expressions is the initial position of 



gas that is at r at time t, and it must be evaluated numerically through the transcendental 



equation (B8). In Figure |4]we compare these analytic predictions to the results of each of the 
magnetized numerical models at t = t en d. The gas density, p, and the non-radial magnetic 
field, Bg, are extracted from the numerical models as azimuthal averages in the midplane of 



the numerical domain where the sine term appearing in equation (B25) is unity. Likewise, 
the radial magnetic field B r is extracted from the numerical models along the x = y 







axis where the cosine term in equation (B23) is unity. The assumption of dynamically weak 
magnetic field is met for r > re in the = 1000 model and we find good agreement between 
the f3 = 1000 model and the analytic prediction at distances not too close to the origin. The 
analytic theory also agrees with the results for stronger fields for r ^ Ar-Q. 



In appendix B. 1 equation (B.2.2 ), we derive an analytic prediction for the total magnetic 



flux that reaches the accretion zone, $ a , under the assumption of dynamically weak magnetic 
fields and neglecting any possible reconnection that occurs near the accretion zone. We have 
assumed that this flux escapes from the accretion zone. Even with reconnection, this method 
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Fig. 4. — Top Left: Azimuthally averaged density in the z = plane. Top Right: Az- 
imuthally averaged radial component of the magnetic field in the z = plane, scaled in- 
versely to the square of r . Bottom Left: Perpendicular component of the magnetic field 
along the x = 0, y = axis, scaled by r and inversely to the density. Each plot shows the 
analytic prediction for the limiting case of weak magnetic field with /3 = 1000. All of the 
plots are taken at time t = t en d- 

accurately tracks the amount of escaping flux, although the time at which the flux escapes 
may be altered by the reconnection. Let $ e sc(^) be the magnetic flux that is inside a radius 
r and that has escaped from the accretion zone. This quantity is well defined only for 
ideal MHD, so that r must be outside the region where magnetic reconnection occurs. At 
large values of r, $ esc ( r ) the total flux released during accretion. As discussed above, 

reconnection occurs in the inner regions of the flow, where it becomes very turbulent. Outside 
this region, the flow is approximately axisymmetric. There we can define ro as the initial 
radius of the gas and magnetic flux, which at time t is located in the midplane at radius 
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r < r . The initial flux inside r is then the sum of the flux inside r(ro) plus the flux that 
has escaped beyond r, 

$o[r (r)] = $(r) + [$ Q - $ CS c(r)], (13) 



where 3>o[?"o( r )] = |B |tt^o- Equation (B8) gives t as a function of r and ro; this can be 



inverted numerically to obtain ro(r, t). We note that equation (13) applies only outside the 
reconnection zone. If we had not assumed that the flux could escape from the accretion zone 
after losing some of its mass, flux would be conserved and both $ a and $ csc would vanish. 

We can use our numerical models to test the predicted value of <3> a and to determine 
the radial distribution of the escaped flux. To do this, we extract $(r) from our numerical 
result at a late time (t = 15£b), and we compare to the analytic result by rewriting equation 



(13) as 

r $q - $esc(r) = $(r) 

$ [r (r)] *„[ro(r)]' 1 ' 

which is the fraction of the escaped flux that is beyond r. In the left panel of Figure 
[5] we show the above expression for the high models. In this case, the assumption of 
dynamically weak magnetic fields used to derive the analytic estimate for r(ro) is well met 
at r ^ re- We expect that $(r) ~ 3>o[ r o( r ')] f° r r ^ r B, and this is confirmed to within 
10% for r > 4r-Q. Given our assumption of a resistive accreting particle, we expect that 
$(r) — > as r — > 0, and Figure [5] confirms this expectation by showing <5$ — > 1 as r — > 0. 
Furthermore, the accumulated flux near the accreting particle shows strong evidence of 
escape for r < 1, consistent with the scale of reconnection-driven tearing modes shown in 



Figure [3] and discussed in section { 3A_ The fact that 5$ is greater than unity at large radii 
is presumably due to the approximation made in determining r (r). In the case of = 100, 
it appears that a significant fraction of the escaped flux ( ^ 20%) has moved outside re- 



In appendix B.l we also predict the radius r$ out to which the magnetic forces associated 



with the accumulated flux strongly affect the flow. The analytic estimate of r$ for high f3 



flow is given by equation (B39 ). In the right panel of Figure |5| we plot this prediction against 
the radius where the median plasma exceeds unity along the perimeter of a control circle in 
the midplane of the high models. At the latest time shown, the prediction agrees with the 
simulation to within about 20% for the = 100 case. Note that at late times, the analytic 
approximation has r$ oc t 1//3 , but it is not known whether the numerical results will continue 
to increase for t > t on( j. It is not entirely clear why the = 1000 results do not agree with 
the approximate model as well as the = 100 results. The model predicts that r$ should 



be very close to (and slightly less than) r$,i, given by equation (B35) for = 1000, whereas 



the simulations show that it is between r$i and r,j, i2 , given by equation (B38). This may be 



associated with the fact that the escaped flux has gone well beyond the sonic point at t en d 



-12- 




Fig. 5. — Left: Radial distribution of escaped flux in the f3 = 100 and j3 = 1000 models at 
time tend- Right: The radial extent of the magnetically dominated region compared to the 
analytic prediction. 

for (3 = 1000 (see the left-hand panel of Figure [5]), so that the conditions are closer to those 
assumed in deriving r$ )2 than for r $1 . 



3.3. Accretion Rate 



Figure [6] shows the rate of accretion onto the central particle as a function of time for 
each of the numerical models. The left plot also includes the result of a purely hydrody- 
namic control model for comparison. As discussed in £ 3.1[ the magnetized models reach 



a statistically steady accretion rate with inflow mediated by reconnection and/or the in- 
terchange instability, whereas the purely hydrodynamic model asymptotically approaches 
the truly steady, spherical Bondi flow. The high frequency modes in Figure [6] have been 
smoothed using a box-car smoothing width of 0.02tB- The red dashed curve shows the 
analytic approximation for the time-dependent accretion rate without magnetic fields from 



appendix B. 1 equation (B15). The analytic estimate is in excellent agreement with the 



purely hydrodynamic numerical model. 

An interesting aspect of the results shown in Figure [6] is that the weak magnetic field 
models (/3 = 100 & = 1000) undergo an initial transient of rapid accretion before settling 
into a steady accretion rate. The reason for this is that enough time must elapse for sufficient 
magnetic flux to accumulate close to the accreting particle for the accretion to the surface of 
the particle to become magnetically dominated, whereas thermal pressure dominates close to 
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Fig. 6. — Accretion rate as a function of time for each of the numerical models compared to 
the analytic prediction for the time-dependent accretion rate for the purely hydrodynamic 
case. In the right plot the time t ss indicated by a gray vertical line when the accretion 
rate is midway between the maximum accretion rate and the final steady state accretion 
rate, representing the characteristic time for the flow to transition from Bondi accretion to 
a magnetically mediated steady state. 



the particle during the initial development of the flow. We can use equation (B8) to estimate 
the time required for the flow to settle into a magnetically mediated steady state accretion 
regime. Specifically, we estimate the time to reach this steady state, t ss , as the time required 



for enough magnetic flux to accumulate inside the thermal sonic radius, r sonic = t-q/2 (Bondi 



1952), so that the average magnetic field within r < r son i C in the midplane corresponds to 
(3 = 1 (i.e., B = (87TP0C 2 ) 1 / 2 for r < r sonic at t = t ss ). Neglecting any flux that has escaped 
beyond r sonic , this then implies 

vrr (r = r sonic , t = t ss ) 2 = vrr 2 onic /3 1/2 (15) 



Solving this for t SB using the transcendental expression for r in equation (B8) determines 
t ss (/3), as shown in Figure [7j The simulations match with this prediction with the (3 = 100 
and /3 = 1000 models transitioning toward the magnetically dominated steady state accretion 
rate at t ~ t ss as shown in Figure [6j 

Figure [8] shows the average accretion rate over the last t& of the simulated time for each 
of the (3 = 10 _1 - (3 = 10 3 models as black circles. The (3 = 10~ 2 model was run only to 
tend = l-5ts and for that case we average over the last £b/2 of the simulated time. The 
vertical bars on each point indicate the standard deviation of the accretion rate over the 
same time interval. It should be noted that these should be interpreted as a measure of the 
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Fig. 7. — An analytic estimate of the time required for enough magnetic flux to accumulate 
inside of the thermal sonic radius for the flow to reach a state of magnetically mediated 
accretion. Black circles indicating the time when the /3 = 10 2 and /3 = 10 3 simulations 
transition from Bondi to magnetically mediated flow are in good agreement with the analytic 
prediction. 

effect of small scale departure from steady accretion flow due to MHD flow instability and 
not as "error bars" in the usual sense of measurement uncertainty. The accretion rate data 
are presented in tabular form as well in Table [2] 

We can obtain a simple analytic model for the accretion flow in the magnetically dom- 



Table 2. Accretion Rates. 
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Note. — Second column: Normalized mean accretion rate for the isothermal equation of 
state models. Third column: Standard deviation of the isothermal accretion rate. 
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Fig. 8. — Average accretion rate as a function of plasma parameter. Error bars show the 
standard deviation in the accretion rate due to interchange and tearing mode unstable flows 



near the accreting particle. The solid line shows equation (17) with the best-fit coefficients 
(3 ch = 5.0 and n = 0.42. 



inated case by assuming that the gas flows in from the Alfven radius tab at the Alfven 
velocity after collapsing vertically from a distance of order the Bondi radius, r B : 

M (x2nr AB -2r AB - Poo v A (xM B (c/v A ) cx M B /3 1/2 ((3 < 1), (16) 

where the second expression follows from equation (B17). We note that |Toropin et al. 



(1999) have shown similar accretion rate dependence with magnetic pressure close to the 
accreteor for the case of the accretion onto a magnetized star. We estimate of the constant 
of proportionality in the above expression that is in rough agreement with our numerical 
results by a more through analytic consideration of the problem in Appendix C.2. Since 



M — > Mb at large /3, a simple relation that captures the limiting behavior in both cases is 



M 



ch 



-i n/2 



-l/n 



+ 1 



(17) 



The solid line in Figure [8] is based on a least-squares fit in log /3 — log M space for the 
parameters /3 c h = 5.0 and n = 0.42. In this notation, /3 c h = 5.0 gives the characteristic value 
of for the transition from the high and low beta limiting cases to occur. 

Sub-grid particle accretion methods have been employed to model protostellar accretion 



in numerical simulations of protostellar cores and clouds by several authors ( 


Bate et al.||1995 


Krumholz et al. 


2004 


Federrath et al. 


2010 


Wang et al. 


2010 


Padoan & Nordlund 


2011) 
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Equation (17) should be of particular utility for extending the sub-grid accretion model for 



embedding Lagrangian sink particles on an Eulerian mesh of Krumholz et al. (2004, 2007) 



and Offner et al. (2009) to the magnetic case for particles moving subsonically through the 



ambient medium. 

It is noteworthy that the qualitative behavior we find at late times is remarkably similar 



to that discovered by Krumholz et al. (2005) for the case of hydrodynamic Bondi accretion 



of a gas with vorticity. The Kelvin circulation theorem for a non-viscous flow is analogous to 
flux- freezing in ideal MHD ( Shu|1992 ), and in the problem of accretion from a vortical fluid, 



the dimensionless vorticity parameter cu* = |V x v|/(c/Vb) defined by Krumholz et al. (2005) 



is analogous to /3~ 1//2 in the present workj^] In both cases, the accretion flow causes a buildup 
of vorticity / flux near the accreting object, which produces regions where the outward 
centrifugal / magnetic force is able to balance gravity and inhibit accretion. For Bondi 
accretion with vorticity, flows with strong vorticity (tu* ^> 1) have steady-state accretion 
rates that scale as roughly / In to* , nearly identical to the /3~ 1//2 scaling we find for the 
strongly magnetized case (/? < 1). For the weak vorticity case (tu* <C 1), the accretion rate 
initially rises to nearly Mb, but then declines as vorticity builds up, reaching an asymptotic 
value < Mb after a transient whose duration is proportional to oo^ 1 . The high /3 cases here 
behave in precisely the same way. 



The only difference we can identify is that, in the vortical case, Krumholz et al. (2005) 



find the accretion rate converges to a value slightly less than Mb in finite time, even in 
the limit — > 0, as long as it is not so small as to place the circularization radius within 
the physical size of the accretor. Here we find that the accretion rate at time t > t ss 
appears to converge to Mb as — > ooj^] The origin of the difference is not entirely clear, 
but one possibility has to do with mechanisms for removing excess vorticity / flux. Both 
can be removed by advection, but magnetic flux can also be rearranged by reconnection, as 
occurs in our simulations. In addition, magnetic buoyancy tends to cause regions of high 



flux to rise away from the accretor. (Similar effects are seen in simulations by Vazquez- 



Semadeni et al. (2011).) In a non-viscous flow, there are no analogous processes capable 



of rearranging the vorticity. In real astrophysical systems, non-ideal MHD and magnetic 
bouyancy effects almost always occur at larger scales than those on which molecular viscosity 
becomes important, and this may lead to a real difference in behavior at late times in the 



1 The —1/2 power arises because the magnetic flux at infinity varies as /3 1 / 2 , while the vorticity at infinity 
scales as cj*. 

2 However, it is not clear from our simulations if the accretion rate would converge to Mb or some lower 
accretion rate in the limit of large /3 at t = oo since then a finite flux could in principle build up near the 
particle. 



weak vorticity / field cases. 



4. Comparison to Adiabatic Models 



It is illustrative to compare our accretion rates to those of earlier works that considered 
the accretion of magnetized gases with a similar field topology but an ideal gas law equation 
of state (7 = 5/3) appropriate for accretion without radiative losses. Pang et al. (2011) 



found that for 1 < (3 < 100, the accretion rate in the adiabatic case depends explicitly on 
the size of the accreting particle with vanishing accretion rate as the particle size — > 0. In 
contrast, for the isothermal case, we find asymptotic convergence toward a finite accretion 
rate with decreasing grid spacing and particle size, even for cases with very strong large scale 



fields (see Appendix |A[) . In the case of adiabatic flow, the results of Pang et al. (2011) and 



Igumenshchev & Narayan (2002) show that mass accumulation in the midplane is limited 



by thermal pressure. In addition, magnetic reconnection leads to thermal pressure-driven 



convective flows that also inhibit mass accumulation in the adiabatic case. The work of Pang 



et al. (2011) has shown that at sufficiently small scale, these effects completely halt accretion. 



Because both of these effects are driven by thermal pressure, neither of them appear in our 
simulations for the isothermal regime. Consequently, radiatively efficient Bondi-type flows 
threaded by large-scale magnetic fields converge to a finite accretion rate in the limit of 
vanishing accreting particle size. 



5. Conclusions 

We have carried out a numerical study of the effect of large-scale magnetic fields in 
an isothermal gas on the rate of accretion onto a resistive point mass — i.e., for the case 
in which only mass, not magnetic flux, accretes onto the point mass. The assumption of 
isothermality is approximately satisfied in regions of star formation, where the cooling time 
of the molecular gas is generally much shorter than the dynamical time for accretion. The 
simulations for this study use simple, very general initial conditions that avoid complications 
arising from boundary conditions by keeping the boundaries far from the accreting object. 
At the same time, our simulations leverage the AMR methodology to retain high spatial 



fidelity close to the accreting object. Contrary to the adiabatic case (Pang et al. 2011), our 
simulations show convergence toward a finite accretion rate as the radiius of the accreting 
object vanishes, regardless of magnetic field strength. We find that magnetic fields reduce 
the Bondi accretion rate in an isothermal medium by about a factor 2 for weak magnetic 
fields (plasma-/? parameter £ 100) at late times, when the magnetic field near the point 
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mass builds up to the point that it can impede accretion. For strong fields (/3 <C 1), the 
accretion rate is reduced by a factor ~ /3 1 / 2 /4. We have developed approximate fitting 
formulae for the accretion rate as a function of /3. The Appendixes give analytic results for 
the time dependent accretion rate of a point mass in the limit of negligible magnetic field 
and for the steady-state accretion rate for the case of a strong magnetic field; both are in 
good agreement with the results of the simulations. 

The authors are grateful of helpful discussions with Eric Agol and Aaron Lee on the topic 
of this paper. Support for this work was provided by: the US Department of Energy at the 
Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344 (AJC and 
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(RIK, CFM, and MRK); the National Science Foundation through grants AST-0807739 
(MRK) and AST-0908553 (RIK and CFM); NSF grant CAREER-0955300 (MRK) and NASA 
through a Spitzer Space Telescope Theoretical Research Program grant (CFM and MRK). 
Support for computer simulations was provided by an LRAC grant from the National Science 
Foundation through TeraGrid resources and the NASA Advanced Supercomputing Division. 
LLNL-JRNL-497719 



A. Numerical Convergence 

The mean steady-state accretion rate at late time is the principle quantity of interest 
from the numerical models presented in this paper. In this section we demonstrate that our 
models provide well-converged estimates for this result. As discussed in §[2j the Alfven-Bondi 
radius tab in our models becomes less resolved as /3 decreases, for a fixed numerical resolution 
scale. Additionally, our numerical models at /3 = 10 3 and (5 = 10 2 use a coarser resolution 
than the lower beta models owing to the computational constraints imposed by the longer 
simulation time required to achieve steady accretion. We therefore focus on demonstrating 
convergence for the set of models that are least resolved in r^B, namely, the model with the 
strongest magnetic field ((3 = 10 -2 ) at resolution Ax = 328 zones/re and the model with the 
strongest magnetic field (/3 = 10 2 ) at the resolution of Ax = 82 zones/re- Figure [9] shows the 
time- dependent accretion rate for each of these models at their native resolution, with the 
resolution and effective sink particle radius coarsened by a factor of 2 and with the resolution 
and effective sink particle radius coarsened by a factor of 4. The convergence properties of 
the instantaneous accretion rate at any particular time is difficult to assess owing to the 
stochastic nature of the accretion rate. However, we can assess the convergence properties 
of the accretion rate averaged over a time interval that is sufficiently long to diminish the 
impact of these stochastic effects. We choose £b/2 for the /3 = 10 -2 model and tb for the 
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Fig. 9.- 
time. 



Convergence properties of the accretion rate of selected models as a function of 



(3 = 10 2 model. 

We find that the late time averaged accretion rates at the resolutions shown in Figure 
M exhibit asymptotic convergence with an implied order of accuracy 



V 



ln(2) 



In 



M 2r - M 



M-M 2x 



(Al) 



that is better than first order accurate. In equation (Al), M is the time-averaged accretion 



rate at the native resolution, M 2x is the time-averaged accretion rate at a resolution that is 
coarsened by a factor of 2 and M± x is the time-averaged accretion rate at a resolution that is 
coarsened by a factor of 4. Given the shock-capturing nature of the RAMSES code, we cannot 
guarantee that better than first order convergence would continue at even higher resolution. 
We therefore estimate the numerical grid convergence error using Richardson-extrapolation 
under the conservative assumption of a first order rate of convergence as 



M-M> 



2.T 



M 



(A2) 



In Table [3] summarizes the convergence properties for each of the models considered in 
this section. We find that the time-averaged accretion rates given by our native resolution 
numerical models are accurate to within 14% of the Richardson-extrapolation estimate of 
the asymptotically converged result. 
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B. Bondi Flow with a Weak Magnetic Field 
B.l. Dynamics 

Here we calculate Bondi flow under the assumption that the gas density is initially 
uniform and then evolves into a steady state. This initial condition corresponds to that 
in our numerical simulations, but would be difficult to realize in practice (for example, an 
approximation to this situation might result when gas flowing supersonically past an object 
is suddenly brought to rest by a strong shock). We assume that the magnetic field is weak 
so that it does not affect the flow. As we shall see below, this approximation breaks down 
sufficiently close to the central mass or at sufficiently late times. The flow is then spherically 
symmetric, and in a steady state the accretion rate is 

M = 47rAr B 2 PocC, (Bl) 

where 

- GM * mo\ 

r B = —^ (B2) 

is the Bondi radius associated with a star of mass M* and A ~ 1.1 for isothermal flow. For 
a steady accretion flow, we then have 

Anr 2 pv = 4nr B 2 pooC. (B3) 

At large radii (r ^> rs), we have p ~ p^ so that 

C T l 

Henceforth, we shall normalize lengths to r&, velocities to c, and times to tb/c; equation 



(B4) then becomes v = r . If we assume that the mass element is initially at rest at Tq, 



then at small radii or at early times, the gas is in free fall, so that 

1 o 1/2 



v = 72 . (B5) 

\r r 



Table 3. Convergence Properties. 



p 


M Ax 


M 2x 


M 


V 


€ 


100 


0.0513 


0.0400 


0.0351 


1.20 


0.14 


0.01 


0.0380 


0.0258 


0.0243 


2.46 


0.062 
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An approximation for the flow everywhere is 



1 (\ 



+ r 2 . 



_ , . , . . (B6) 

v y/2 \r r / 

It should be noted that, although we used the approximation of steady flow to estimate the 



-1/2 



velocity at large radii, equation (B6) for the velocity is time dependent: r is a function of 



both r and t, so dv/dt ^ 0. In equation (B15) below we shall give the time-dependent result 



for Bondi flow that occurs in an initially stationary medium. 

How long does it take a particle to reach a point r when it starts at ro? Integration of 



equation (B6) gives 
t = 



r ° dr 
v 



_3/2 



[x(l — x)] 1 ^ 2 + arctan 



1 — x 



1/2' 



r%{l-x 3 ), 



X/2 I ' \ -r 

where x = t/tq. The time at which the gas is accreted at the origin (x = 0) is 



t a 



2 3/2 



(B7) 
(B8) 

(B9) 



Note that this result is approximate, since it depends on the harmonic mean approximation in 



equation (B6). We have found better agreement with the numerical results if we approximate 



t a as the root mean square of the two terms in equation (B9): 

1/2 



t n 



n \,3 , 1 6 



9 



(BIO) 



The solution of this equation shows that gas accreting at time t originated from a radius r 0a 
given by 

^ (Bll) 



0a 



where 



16 J l + (l + r 2 )V2 : 
/ 16 



r S ^l* = 0.640t 



For late times (t ^> 1), this reduces to 



r 0a -> (St) 



1/3 



(B12) 



(B13) 



The accretion rate onto the origin is 

M = 47rAr 2 aPc 



dr, 



dt 



Oa 2 

— x r B c, 



(B14) 
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where the final factor gives M the correct dimensions. Evaluating the time derivative from 



equation (Bll), we obtain an approximation for the time-dependent accretion rate, 

M(t) ~ 4vrr B 2 pooC 



1 + r 



2^/2' 



(B15) 



Thus, at early times the accretion rate increases linearly with time, whereas at late times it 
approaches the steady state value given in equation (Bl ) (although here we have set A = 1). 



Consider now the particular case of steady flow. Since the initial location of a mass 
element, t*o, depends on both r and t, the steady flow approximation is valid only if the 



1 term in equation (B6) is negligible. This is true for r < r or for sufficiently large r$ 



provided r is not too close to r . As a check on the accuracy of equation (B6) in this case 
(i.e., when r$ 1 is negligible), note that the actual sonic point is at re/2 (Shapiro & Teukolsky 



1983), whereas equation ( |B6[ ) gives 0.65re; the approximation is thus accurate to within 



about 30%. Equation (B3) gives the density for a steady flow, which requires that the r, 



term in equation (B6) be negligible: 

P 1 , 



(steady flow). 



(B16) 



B.2. The Magnetic Field 

When gas accretes onto the central object, both its mass and its pressure are removed 
from the ambient medium. In the case of the magnetic field, we assume that the flux is not 
accreted by the central star. As a result, the flux associated with the accreted matter, $ a , 
builds up and distorts the flow close to the central object. When a flux tube loses mass, it 
becomes buoyant and drives an interchange instability. However, gas continues to accrete 
along this flux tube so it may eventually fall back to the center. We therefore expect that 
the innermost region will become turbulent. We begin with a discussion of the magnetic 
field in the absence of the effects of the accretion flux, and then estimate its effect at the 
end. 



B.2.1. The Field in Smooth Inflow 

Just as the gravitational force due to the star becomes important at radii less than the 
Bondi radius, r B , in the hydrodynamic case, so we expect it to become important at radii 
less than the Alfven-Bondi radius, 

GM» 4ttGM, Poo 

TAB = —2~ = 7^2 , (B17) 
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,15 M*/M Q 



= 3.32 x 10 15 - , 7 - 1v> cm (B18) 

(u A /2 km s" 1 ) 2 V ; 



in the MHD case. The ratio of the Alfven-Bondi radius to the standard Bondi radius is 

TAB C 2 1 n 



2 



2 P ' 



(B19) 



where (3 = 8n p^c 2 / B 2 is the plasma (3. Our assumption that the field is weak implies (3^1. 
There is an important relation between tab and the magnetic critical mass 

M. = ^ (B20) 

which also determines the relative importance of self-gravity and magnetic fields: 

r AB 4ttGM, Poo 3fM MA 

^ = ^r = A~w)^ (B21) 

where M = 47rp oo r ! /3. The magnetic field is dominant for r > tab- in the purely gaseous 
case, the mass is subcritical for M < M$; in the Bondi case, we see that the gas mass M is 
replaced by the geometric mean of the gas mass and the stellar mass (ignoring the factor |). 
Shu, Li, & Allen (2004) obtained a similar result for the case in which the gas is in a disk; 
they showed that it was possible for the field to be so strong that it could "levitate" the gas 
above a star in the process of formation. Note that the fact that it is the geometric mean 
mass that determines whether the gas is sub- or super-critical has an important consequence: 
in the purely gaseous case, a sufficiently large uniform cloud is always supercritical, since 
Mo oc rj] and $ oc r 2 . However, in the Bondi case, the opposite occurs: a sufficiently large 
cloud is always subcritical, since now (MoM*) 1 / 2 oc r^ 2 increases more slowly than $. 

We assume that the field is initially uniform, so that B^q = 0; for spherical inflow, B^ 
will remain zero. For a spherical inflow, the radial flux through any surface r 2 dQ remains 
constant, so that 

r 2 B r dVL = r 2 B r0 dtt, (B22) 



which implies 

B r = I )' r o 

r / V r 



ey^B^On 2 . (B23) 



To evaluate Bg, consider a spherical shell of thickness dr and radius r. The flux in the 
shell at 9 is proportional to Bgrdr. The mass in the shell is 47rr 2 pdr. Since each of these 
remains constant in the inflow, we have 

rB e dr oc pr 2 dr, (B24) 
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which implies 

B e = B eo ( = -B sin 9 ( , (B25) 
where the sign corresponds to the case in which the initial field is Bo = Bqz. 

How does the magnetic force compare with the gravitational one? First, we note that 
the radial field by itself exerts no force; we therefore consider the pressure exerted by Bg 
and the tension force. We consider times late enough so that r t ~ (3t) 1//3 and thus that r$ 
is approximately independent of r. For the pressure force, the relative importance of the 
magnetic field and gravity in the midplane (9 = vr/2) can be assessed from the ratio 

v\ Bgr ( p \ r 3 
4 ~ Att P GM ~ V/W r AB r 2 ' 

At large radii, we have p ~ p^ initially (r ~ tq) the magnetic field dominates for r > tab, 
as expected. At small radii, p/poo oc r~ 3 / 2 so that magnetic effects ex v\jv\ oc r 3//2 become 
negligible. 

Next consider the tension in the radial direction, 
The ratio of this force in the midplane to the gravitational force is 

tension (B VB) r T ^ ^ ^ ^ 

(B29) 



F g AnGMp/r 2 r AB V Po^ 3 , 

r 



tab 



r 3 / 1 
1 =-11 



V^r 3 / 2 



where the last expression applies to steady flows. Provided r ^ 1, the density dependent 
term becomes negligible for r -C ro, so that in this case the force ratio becomes 



^tension _ '"O (B30) 



F g r AB 



Since r$ ~ (St) 1 / 3 at late times (eq. B13), it follows that the tension force will eventually 
dominate and render the accretion anisotropic at 



(B3i) 

where we have explicitly included the factors of r B and c. This is to be expected, since as 
noted above a sufficiently large cloud is subcritical. 
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B.2.2. Effects of the Accretion Flux 

The accretion flux, <3> a , is the magnetic flux associated with the mass that has accreted 
onto the central mass. We expect this flux to be buoyant and to therefore lead to turbulence. 
Here we estimate the size of region affected by the accretion flux. 

At a time t, the accretion flux is the flux inside the initial radius r^ a given in equation 



(BID 



<f> /n 2\2/3 4/3 

" 4=(^r) Z J z^m , (B32) 



nr B 2 B 0a \16j [1 + (1 + r 2)i/2]V3 ; 

We estimate the radius, r$, out to which this flux extends by assuming that the field asso- 
ciated with $ a is uniform, and that the flow at r<j, is steady. The latter assumption requires 



that r$ be small compared to the starting radius, r , since as discussed below equation (B6), 
r (r,t) introduces time dependent effects. We consider two limiting cases: (1) r$ <C 1, where 
the accretion flux interacts with supersonic inflow and (2) r$ > 1, where the accretion flux 
interacts with the pressure in the ambient medium. 

Case 1: Supersonic inflow (early and intermediate times): We estimate r^i, the value 
of r$ in this case, by determining where the pressure due to the accretion field balances the 
ram pressure of the accreting gas. Since we are assuming that r $1 < 1 and r$ i <C r , 



equations (B6) and (B16) imply 



B\ 2 a/2PooC 2 /T3QQ\ 

8^ =pV = (B33) 



Flux conservation implies B a r\ x = B rl a , so that 

„8/3 



= ^fe^ (B34) 



21/3^2/3- 

1 /9tt 2 \ 8/9 r 16 / 9 



2V3/32/3 ^ 16 J [1 + (1 + r 2)l/2]^ 



(B35) 



This expression is valid for both r < 1 (early times) and r > 1 (intermediate times). At late 
times, the flow is dominated by thermal pressure. 

Case 2: Pressure- confined flow (late times): In this case the magnetic pressure associ- 
ated with the accretion flux balances the thermal pressure of the ambient medium, 

g = ^ * (B36) 
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Flux conservation then implies 



r *,2 = (B37) 

1 ^9tt 2 \ 1/3 r 2 / 3 , 
^173 - (B38) 



[l + (l +r 2)l/2]V3- 

In order to obtain an approximation valid at all times, we write 

1 1 

r* ~~ r* i 




(B39) 



Note that r$ is less than either r $1 or , r$ 2 corresponding to the fact that in this simple 
model the pressure due to the escaped flux has to balance both the thermal pressure and the 
ram pressure. Since r% 1 /r% 2 exceeds unity only at late times, this can be approximated as 

3.6r 16 / 9 1 ,„ , 

r$ ~ • (B40) 

^ 2 / 3 [l + (l + r 2 )V2] 8 / 9 (l + 4.0/3- 5 /6 r io/9)V2 • 

At early times, r$ oc r 16 / 9 ; at intermediate times (1 <Cr< 0.3/3 3 / 4 ), r$ oc r 8 / 9 ; and at late 
times r,j> oc t 1 / 3 . 



C. Magnetic Bondi Flow in a Strong Magnetic Field 

C.l. Initial Transient 

A striking feature of Figure 2 for strong fields is that the flow is isotropic beyond some 
radius, but then predominantly aligned along the axis inside that, until the flow is very close 
to the center. This makes sense, since initially the field is straight and therefore exerts no 
force; thus, at sufficiently early times, the flow for a strong field is almost identical to that for 
no field. We focus on the region inside re, where we neglect pressure forces. Let r = r — 5, 
where 5 r since we are considering early times. Then equation (B5) implies 

d5 [2r B 5\ 1/2 (2r B( 5)V 2 

v = — = c ~c , (CI) 

dt \ rr J r 

where we have written the equation in dimensional form. Integration gives 
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The ratio of the tension force to the gravitational force at early times is given by equation 
(B28) with p = p . For small 8, this is 

# = (C3) 
F„ r AB ^ ' 

The magnetic field will begin deflecting the flow from a radial trajectory to an axial one 
when this ratio is of order unity, which occurs at 

We have found that the growth of the region deflected from a radial trajectory in our nu- 
merical simulations with /3 — 0.1 and (3 = 0.01 follow this functional form very well but that 
the deflection from spherical flow occurs somewhat later than predicted. We extract a good 
empirical fit to the low-/? simulations with 

r /2.0\ 1/2 t 



C.2. Magnetic Bondi Flow in a Strong Magnetic Field ((3 £ 0.1) at Late Times 

For a very strong field, the gas will attempt to settle into vertical hydrostatic equilibrium, 

P = Pooe- m * /kT = Pooe rB/r (C6) 

where m is the mass per particle and = —GM*/r is the gravitational potential. Henceforth, 
we shall normalize all lengths to the Bondi radius, as in the previous section. Outside the 
Bondi radius, this expression gives only a modest increase in density, but for small radii the 
increase can be very large-so large that it takes a long time to reach equilibrium. Let w be 
the cylindrical radius, so that r = (w 2 + z 2 ) 1 ^ 2 , where z is the height above the disk. The 
density at the midplane (r = zu) is then 

Po = Pooe 1/ro . (C7) 

For small radii, nj<Cl, the density distribution near the midplane is approximately gaussian, 

p^poe"^ 2 , (C8) 

where p is the midplane density and the scale height is 

h = v/2w 3 / 2 . (C9) 



-28- 



In equilibrium, the total surface density of the gas near the midplane is then 

S eq ~ 2p h = Poo r B (2w) 3/ V /ro , (CIO) 



where we have used equation (C7) to eliminate p . 



When do magnetic forces balance gravity? For a thin disk, magnetic tension dominates 



magnetic pressure (Shu & Li 1997). For an axisymmetric field, the net radial tension is 



1 



(B • V B t 



1 



4/T x ' 47itb * dz 

Integrating through the disk, we find that the forces balance when 

—B z (2B m ) = 

where B m is measured just above the disk. 



(Cll) 



(C12) 



To obtain an accurate solution beyond this point, we would have to solve for the struc- 
ture of the field. This is a challenging problem even when the system is in equilibrium. Here, 
however, we are assuming that the system is in equilibrium outside some critical radius, w CT , 
but that there is an unknown accretion flow inside that radius. We therefore content our- 
selves with attemping to infer the scaling for the solution. We assume that B z in the disk 
is proportional to the ambient field, B^, and that the radial component of the field, _B ro , is 



proportional to B z . Equation (C12) then implies that 



w 



(C13) 



For a given location in the disk, gas will accrete along the field lines until the surface density 
reaches this value. The field is unable to support more gas than this, so this value represents 
an upper limit on E; any additional gas will accrete onto the central star. However, we 



have determined another maximum value for the surface density in equation (CIO), which 



is the value the surface density has in hydrostatic equilibrium. Equating these two surface 
densities determines the critical radius, zu CT : The gas can be supported by the field outside 



zu CI , but inside w cr gas that exceeds the surface density in equation (C13) must fall onto the 



central star. Equations (CIO) and (C13) imply that this critical radius satisfies 



P- 



(C14) 



A good approximation for the solution of this equation for (5 0.15, corresponding to 



0.6, is 



ln/3" 1 -0.5 In In 



(P < 0.15). 



(C15) 
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In the regime of greatest interest, 10 3 < f3 < 0.15, the solution can be approximated by the 
simpler form 

zu CI ~ 0.85/3 1/4 (10~ 3 < (3 < 0.15). (C16) 

The accuracy of this solution in the prescribed range is about 10%, which is much better 
than the accuracy of the underlying equation. 

We are now in a position to estimate the accretion rate onto the central star. We 
assume that the accretion flow is primarily along the field lines, and that it is initiated by a 
rarefaction wave propagating at the sound speed, c. After an initial phase during which the 
surface density just inside zu cr becomes large enough that it distorts the field so much that 
it can accrete, the accretion rate on both sides of the disk becomes 



M~2(7rrB 2 ^ iCr ) Poo c, (C17) 

where cCoo.cr is the cylindrical radius of the critical field lines far from the star. If we assume 
that zu CI oc t£7oo,er> then in the range 10~ 3 £ (3 £ 0.1 we have M oc w 2 cx oc /3 1//2 , and we can 
write 

M = 47rA low/3 r B 2 pooC/3 1/2 , (C18) 

where Ai ow /3 is a numerical constant. Note that the /3 1//2 scaling is the same as that implied 
by the crude argument in the text. Were we to assume that ztfoo.cr — w a and that equation 



(C16) were accurate, then Ai ow ^ would equal 0.36. This estimate is within a factor 1.6 of the 
numerical results. Setting Ai OW/ 3 = 0.24 gives an accretion rate that agrees with the results 
of the simulations for (3 = 0.1, 0.01 to within 8%. 
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